library("plyr")
library("ggplot2")
library("fUnitRoots")
library(tseries)
library(forecast)
#读取数据
data<-week_de("product_1")
data<-week_de("product_2")
data<-week_de("product_3")
#筛选数据
tmp<-data$p$all-data$r$all
tmp<-tmp[data$pr$full=="h"&data$pr$median=="h"&data$pr$zero=="h",]
#tmp<-tmp[data$pr$full=="l",]
sum(tmp)
nrow(tmp)

#日刊预测
ts<-lgts(tmp,n=8,start=50,end=102,k=4)
start=63
#经测试51周为周期效果最佳(增量时51，减量时50)
end=start+51
res<-myaim(data=tmp,id=6,start,end,k=4)
#debug(myaim)
res[[1]]
res[[2]]
rdata2<-res[[3]]
summary(rdata2)
mod<-rdata2$model
mod$phi
mod$theta
mod$Delta
debug(myaim)
#期刊预测
data<-tmp[1,]
data<-data[data!=0]
start = 1
end = 24
data_t<-data[((start-1)*k+1):(end*k)]
#下面是周期测试函数
i=16
for(i in 2:30){
wts<-ts(data[data!=0],start = 1,end=c(floor(90/i),i),frequency = i)
timed_t<-decompose(wts)
print(max(timed_t$seasonal)-min(timed_t$seasonal))
}
plot(timed_t)
data<-tmp[4,]
data<-data[data!=0]
wts<-ts(data,start = 1,end=c(floor(100/13),13),frequency = 13)
wts[wts>(3*sd(wts)+mean(wts))]<-mean(wts)
pp<-list()
pp[[1]]<-decompose(wts)
plot(pp[[1]])
lg<-plts(ts=pp,lag.max=40,lgval=0.15,dif=4,sen=FALSE,diff=0)
lg[[2]]
lg<-plts(ts=pp,lag.max=40,lgval=0.15,dif=13,sen=TRUE,diff=0)
lg[[2]]
debug(plts)
undebug(plts)
debug(myaim_wk)

res<-myaim_wk(tmp,id=4,start=3,end=6,k=13)
res[[1]]
#debug(myaim_wk)
res[[2]]
summary(res[[1]])
rdata<-res[[1]]
rdata$call
#debug(predict(res[[1]],n.ahead=rd))
predict(res[[1]],n.ahead=rd)
rdata$model$phi
rdata$model$theta
rdata$model$Delta
rdata$model$Z
rdata$model$a
rdata$model$P
rdata$model$T
rdata$model$V
rdata$model$Pn
rdata$residuals
,D=0,start.P = 2,start.Q = 1,max.p=12,max.q=2,trace=Tmax.p=12,max.q=2,start.p=12,start.q=0,
wts<-ts(tmp[1,],start = 8,end=c(14,8),frequency = 8)
arima1<-auto.arima(wts, ic = c("aic"),D=1,trace=T,seasonal = TRUE)
attach(rdata$model)
rdata<-res[[1]]
rdata<-rdata$model$T%*%rdata$model$a+rdata$model$V%*%rdata$model$a
res[[2]]
(res[[3]][length(res[[3]])]-res[[3]][length(res[[3]])-13])*0.22662

forecast<-forecast.Arima(ARIMA,h=13,level=c(99.5))
debug()
debug(forecast.Arima)
plot.forecast(forecast)

mmean<-function(data=tmp,id=1,n=13,start=12,end=23,k=4){
  names<-rownames(data)[id]
  data<-data[id,data[id,]!=0]
  data_t<-data[((start-1)*k+1):(end*k)]
  wts <- ts(data_t,start = start,end=c(end,k),frequency = k)
  wts[wts>(3*sd(wts)+mean(wts))]<-mean(wts)
  ts <- decompose(wts)
  #plot(ts)
  fit<-ts$x
  last<-length(fit)
  for(i in 1:12){
    fit<-c(fit,mean(fit[(last-n+1):last]))
    last<-length(fit)
  }
  fit
}

mean((x1-xn)/n)


fit<-mmean(data=tmp,id=1,n=13,start=12,end=23,k=4)
plot(fit,type="l")




